function [sgreduced, ...
          id_Ugridreduced, s_id_gridreduced, sU_id_gridreduced,...
          P_col,...
          Tcoord_zeros,...
          Ugridreduced, n_U_sets, indices_pairs_extended, Tcoord, C, C_sign,...
          cr_obsequiv, cr_one, cr_zeros,fill_obsequiv, fill_one, fill_zeros, numeqconstraints, RHS_one, RHS_zeros,RHS_obsequiv,...
          equalorder, n_vertices,infty_pos,...
          cr_symmetry, fill_symmetry, RHS_symmetry, S1, S2, S3, ind,...
          sU_id_gridreduced_2]=useful_anyparam(u0grid, u1grid, u2grid, u3grid, u4grid, sg,random_indices_pairs_old, random_indices_pairs_new,PYX)           
%% U sets
n_U_sets=4; %number of U-sets

Ugrid{1}=[u0grid-u1grid   u0grid-u2grid  u0grid-u3grid  u0grid-u4grid u1grid-u0grid...
          Inf*ones(sg,1) -Inf*ones(sg,1) ...
          -(u0grid-u1grid)    -(u0grid-u2grid)   -(u0grid-u3grid)   -(u0grid-u4grid)  -(u1grid-u0grid)...
          zeros(sg,1)];

Ugrid{2}=[u2grid-u1grid   u1grid-u2grid  u1grid-u3grid u1grid-u4grid u2grid-u0grid ...
          Inf*ones(sg,1) -Inf*ones(sg,1) ...
           -(u2grid-u1grid)    -(u1grid-u2grid)   -(u1grid-u3grid)  -(u1grid-u4grid)  -(u2grid-u0grid) ...
          zeros(sg,1)];
      
Ugrid{3}=[u3grid-u1grid   u3grid-u2grid  u2grid-u3grid u2grid-u4grid u3grid-u0grid ...
          Inf*ones(sg,1) -Inf*ones(sg,1) ...
           -(u3grid-u1grid)    -(u3grid-u2grid)   -(u2grid-u3grid)  -(u2grid-u4grid)  -(u3grid-u0grid) ...
          zeros(sg,1)];    

Ugrid{4}=[u4grid-u1grid   u4grid-u2grid  u4grid-u3grid u3grid-u4grid u4grid-u0grid ...
          Inf*ones(sg,1) -Inf*ones(sg,1) ...
           -(u4grid-u1grid)    -(u4grid-u2grid)   -(u4grid-u3grid)  -(u3grid-u4grid)  -(u4grid-u0grid) ...
          zeros(sg,1)];  
%% Reduce number of grid points to check by applying our partitioning arguments
%equalorder (sgx1) assigns the same index to parameters vectors such that Ugrid{1}, Ugrid{2}, ..., Ugrid{n_U_sets} have the same order

d=cell(n_U_sets,1);
index=cell(n_U_sets,1);

for j=1:n_U_sets 
    [sortedUjgrid,index{j}] = sort(Ugrid{j},2); 
    %sortedUjgrid gives Ugrid{j} with each row sorted from smallest to largest; it is sgxn_columns_Ujgrid
    %for each row i of sortedUjgrid and for k=1,2,...,n_columns_Ujgrid, index{j}(i,k) gives the column position in Ugrid{j}(i,:) of sortedUjgrid(i,k); it is sgxn_columns_Ujgrid
    d{j} = diff(sortedUjgrid,[],2) == 0; %sg x (n_columns_Ujgrid-1)
    %Suppose that Ugrid{j} has 3 columns. Then, for each row i of sortedUjgrid, d{j} will be: 
    %[1 1] if sortedUgridj(i,2)-sortedUgridj(i,1)=0 and sortedUgridj(i,3)-sortedUgridj(i,2)=0
    %[1 0] if sortedUgridj(i,2)-sortedUgridj(i,1)=0 and sortedUgridj(i,3)-sortedUgridj(i,2)~=0
    %[0 1] if sortedUgridj(i,2)-sortedUgridj(i,1)~=0 and sortedUgridj(i,3)-sortedUgridj(i,2)=0
    %[0 0] if sortedUgridj(i,2)-sortedUgridj(i,1)~=0 and sortedUgridj(i,3)-sortedUgridj(i,2)~=0
end

[~,~,equalorder] = unique([[index{:}] [d{:}]],'rows', 'stable'); %sgx1 

%Reduced grid of parameter vectors and U sets
[~,index_final,~]=unique(equalorder, 'stable'); %one row from equalorder per ordering group together with its index
Ugridreduced=cell(1,n_U_sets);
for j=1:n_U_sets 
    Ugridreduced{j}=Ugrid{j}(index_final,:); %sgreduced x n_columns_Ujgrid
end
sgreduced=size(Ugridreduced{j},1); %size of the reduced grid of parameter values


%% Construct indices to keep track of equal elements inside Ugridreduced{j} for j=1,...,n_U_sets
id_Ugridreduced=cell(n_U_sets,1);
for j=1:n_U_sets 
    id_Ugridreduced{j} = NaN(size(Ugridreduced{j})); % %sgreduced x n_columns_Ujgrid
    for k = 1:sgreduced
        [~, ~, id_Ugridreduced{j}(k,:)] = unique(Ugridreduced{j}(k,:), 'stable'); 
    end
end

s_id_gridreduced=zeros(sgreduced, n_U_sets); %useful to set number of unknowns of the linear programming
for j=1:n_U_sets
    s_id_gridreduced(:,j)=max(id_Ugridreduced{j},[],2); %number of distinct elements of Ugridreduced{j}
end
sU_id_gridreduced=prod(s_id_gridreduced,2); %sgreducedx1 reporting the cardinality of the Cartesian product of the U sets
                                            %sU_id_gridreduced(g) is the number of unknowns of the LP for each g-th parameter value
s_id_gridreduced_2=zeros(1, n_U_sets); 
for j=1:n_U_sets
    s_id_gridreduced_2(:,j)=size(id_Ugridreduced{j},2); %number of elements of Ugridreduced{j} (count as separate even equal elements)
end
sU_id_gridreduced_2=prod(s_id_gridreduced_2); %cardinality of the Cartesian product of the U sets 

%% Invariant components of the linear programming problem
infty_pos=6; %position of Infty in U-sets
mininfty_pos=7; %position of -Infty in U-sets
zero_pos=13;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Equalities %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Observational equivalence
% Create matrix P 15x4 such that
% the first row of P is [1 1 1 1]
% from the 2nd to the 5th of P we have 1 only once per row
% from the 6th to the 11th of P we have 1 twice per row
% from the 12th to the 15th of P we have 1 three times per row
% The total number of rows of P comes from the number of elements of each observational equivalence constraint
% The total number of columns of P is equal to n_U_sets 
P_temp = dec2bin(1:2^n_U_sets-2)-'0'; 
[~, ind] = sortrows([sum(P_temp,2) -P_temp]); 
P = [ones(1,n_U_sets); P_temp(ind,:)]; 
length_P=size(P,1);
% Coordinate columns
P_col=cell(n_U_sets+1); 
for j=1:n_U_sets+1
    P_temp=P;
    P_temp(P==0)=infty_pos; %infty
    P_temp(P==1)=j;
    P_col{j}=P_temp;
end

% Coordinate rows
cr_obsequiv=cell(n_U_sets+1); 
for j=1:n_U_sets+1
    cr_obsequiv{j}=j*ones(1,length_P);
end
cr_obsequiv=cell2mat(cr_obsequiv).';
cr_obsequiv=cr_obsequiv(:).';

% Fill
fill_obsequiv=cell(n_U_sets+1); 
for j=1:n_U_sets+1
    fill_obsequiv{j}=[1 -ones(1,4) ones(1,6) -ones(1,4)]; %specific (opposite to notes)
end
fill_obsequiv=cell2mat(fill_obsequiv).';
fill_obsequiv=fill_obsequiv(:).';

RHS_obsequiv=(-1+PYX).';
%% Zeros
% Create matrix Tcoord of dimension sU_id_gridreduced_2xn_U_sets reporting the column coordinates of all the possible n_U_sets-tuples from the U-sets
delta=6; %extra elements of U-sets with respect to Restrictions1
%vectors = arrayfun(@(x) {1:x}, s_id_gridreduced_2-delta); %Tcoord as for Restrictions1
%n = numel(vectors);
%Tcoord_temp = cell(1,n);
%[Tcoord_temp{:}] = ndgrid(vectors{:}); 
%Tcoord_temp = cat(n+1, Tcoord_temp{:});
%Tcoord_old = reshape(Tcoord_temp,[],n); 

vectors_augmented = arrayfun(@(x) {1:x}, s_id_gridreduced_2);  %Tcoord as for Restrictions2
n = numel(vectors_augmented);
Tcoord_temp = cell(1,n);
[Tcoord_temp{:}] = ndgrid(vectors_augmented{:}); 
Tcoord_temp = cat(n+1, Tcoord_temp{:});
Tcoord_augmented = reshape(Tcoord_temp,[],n); 

aug_size = ones(1,n)*delta; % augment size of each vector
vectors_aux = arrayfun(@(a) {[false(1,s_id_gridreduced_2(1)-delta) true(1, a)]}, aug_size);
T_aux = cell(1,n);
[T_aux{:}] = ndgrid(vectors_aux{:}); 
T_aux = cat(n+1, T_aux{:});
T_aux = reshape(T_aux,[],n); 
[~, ind] = sortrows(T_aux, n:-1:1); % indices of stably sorting the rows.
% Most significant column is rightmost, as per your code
Tcoord = Tcoord_augmented(ind, :);

Tcoord_zeros=Tcoord(any(Tcoord==mininfty_pos,2),:); %rows of Tcoord with at least one mininfty_pos (-Inf) 
n_zeros=size(Tcoord_zeros,1);
cr_zeros=cr_obsequiv(end)+1:1:cr_obsequiv(end)+n_zeros;
fill_zeros=ones(1,n_zeros);
RHS_zeros=zeros(n_zeros,1); 

%% Integrates to one
cr_one=cr_zeros(end)+1;
fill_one=1;
RHS_one=1; 

%% Symmetric marginals (other nonparametric distributional restrictions and zero volume boxes already imposed above)
% Create matrix S1 of dimension ((n_U_sets+1)*n_U_sets) x (n_U_sets) which lists the positions of the elements entering the first part of the symmetry sum
S_temp=cell(n_U_sets+1,1);
for i=1:n_U_sets
    S_temp{i}=ones(n_U_sets+1,n_U_sets)*infty_pos;
    S_temp{i}(:,i)=(1:1:n_U_sets+1).';
end
S1=vertcat( S_temp{:});

% Create matrix S2 of dimension ((n_U_sets+1)*n_U_sets) x (n_U_sets) which lists the positions of the elements entering the second part of the symmetry sum
S_temp=cell(n_U_sets+1,1);
for i=1:n_U_sets
    S_temp{i}=ones(n_U_sets+1,n_U_sets)*infty_pos;
    S_temp{i}(:,i)=(mininfty_pos+1:1:mininfty_pos+n_U_sets+1).';
end
S2=vertcat( S_temp{:});

% Create matrix S3 of dimension n_U_sets x n_U_sets which lists the positions of the elements entering the zero symmetry constraint
S3=eye(n_U_sets)*zero_pos;
S3(S3==0)=infty_pos;

cr_symmetry=cr_one(end)+([kron(1:1:(n_U_sets+1)*n_U_sets, ones(1,2)) (n_U_sets+1)*n_U_sets+1:1:(n_U_sets+1)*n_U_sets+n_U_sets]);
fill_symmetry=ones(1,size(cr_symmetry,2));
RHS_symmetry=[ones((n_U_sets+1)*n_U_sets,1); 1/2*ones(n_U_sets,1)];

numeqconstraints=(n_U_sets+1)+1+n_zeros+((n_U_sets+1)*n_U_sets+n_U_sets);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Inequalities %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Increasingness 
% Create matrix indices_pairs_extended of dimension %[sU_id_gridreduced_2*(sU_id_gridreduced_2-1)/2]x[2] reporting the 
% row indices of all possible unordered pairs of rows from the matrix obtained by taking the Cartesian product of the U sets
%indices_pairs_extended=pairIndices(sU_id_gridreduced_2);
%Continuing the example above when constructing Tcoord
%indices_pairs_extended=[1 2; ...; 1 27; 2 3; ...; 2 27; ... 26 27]

old_length=7;
[C1, C2] = myind2ind (random_indices_pairs_old, old_length^n_U_sets);
indices_pairs_extended_old=[C1 C2];
[C1, C2] = myind2ind (random_indices_pairs_new, sU_id_gridreduced_2);
indices_pairs_extended_new=[C1 C2];
indices_pairs_extended=[indices_pairs_extended_old; indices_pairs_extended_new];

% Create matrix C of dimension 2^(n_U_sets)xn_U_sets from A=[1:1:n_U_sets] and B=[n_U_sets+1:1:2*n_U_sets] such that for each row
% first element is A(1) or B(1)
% second element is A(2) or B(2)
% ...
% last element is  A(end) or B(end)
A = 1:1:n_U_sets;
B = n_U_sets+1:1:2*n_U_sets;
t = dec2bin(0:2^n_U_sets-1)-'0';
[~, ind_sort] = sortrows([sum(t,2) -t]);
t = t(ind_sort, :);
AB = [A B];
ind_AB = t*n_U_sets + (1:n_U_sets); 
C = AB(ind_AB);

% Create matrix C_sign of dimension 1x2^(n_U_sets) where C_sign(i)=1 if C(i.:) contains an
% even number of elements from A, -1 if C(i.:) contains an odd number of elements from A
C_sign_temp=mod(sum(C<=n_U_sets,2),2);
C_sign=C_sign_temp;
C_sign(C_sign_temp~=0)=-1; %if odd
C_sign(C_sign_temp==0)=1; %if even
C_sign=C_sign.'; %1x4

%% Other
n_vertices=2^n_U_sets;
end